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Abstract 

We analyze the stability of a cylindrical Couette flow under the impo- 
sition of a weak axial flow in case of a very short cylinder with a narrow 
annulus gap. 

We consider an incompressible viscous fluid which is contained in the 
narrow gap between two concentric short cylinders, where the inner cylin- 
der rotates with constant angular velocity. The caps of the cylinders have 
narrow tubes conically tapering to super narrow slits which allow for an 
axial flow along the surface of the inner cylinder. 

The approximated solution for the Couette flow for short cylinders was 
found and used for the stability analysis instead of the exact but bulky 
solution. The sensitivity of the Couette flow to general small perturbations 
and to the weak axial flow was studied. We demonstrate that perturbations 
coming from the axial flow cause the propagation of dispersive waves in 
the Taylor-Couette flow. 

The coexistence of a rotation and of an axial flow requires to study 
in addition to the energy and the angular momentum also the helicity of 
the flow. The approximated form for the helicity formula in case of short 
cylinders was derived. 

We found that the axial flow stabilizes the Taylor - Couette flow. The 
supercritical flow includes a rich variety of vortical structures including a 
symmetric pair of Taylor vortices, an anomalous single vortex and quasi 
periodic oscillating vortices. Pattern formation was studied at large for 
rated ranges of azimuthal and axial Reynolds numbers. A region where 
three branches of different states occur was localized. Numerical simula- 
tions in 3D and in axisymmetrical case of the model flow are presented, 
which illustrate the instabilities analyzed. 
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1 Introduction 

Since the famous experiments of Taylor £Q in long cylinders many experimental 
and theoretical investigations of the interesting phenomena of arising and evolu- 
tion of Taylor vortices were done. The main part of these studies was devoted 
to cases of long cylinders. It was done under the assumption that a sufficiently 
long cylinder and periodic boundary conditions will emulate the infinitely long 
cylindrical region well. In addition such regions are convenient for theoretical 
investigations. The experimental investigations of short cylinders started much 
later and gave surprising results (see, for example, Benjamin and Mullin [2J), 
which indicated that the zone and the magnitude of the influence of the caps of 
the cylinders is much bigger and takes effect on the type of motion in the whole 
region. It was assumed that close to the caps of the cylinders the boundary layer 
is directed from the outer to the inner cylinder as a result of the existence of 
an Ekmans boundary layer. However, the experiments gave a more complicated 
picture of the motions. There exists also an atypical layer in opposite direction, 
i.e. from the inner to the outer cylinder connected withe the existence of the 
'anomalous mode'. The first mathematical model which took the influence of 
the cylinder caps into account was represented in the work of Schaeffer This 
model is not very consistent and needs improvements. The role of the anomalous 
mode was studied later numerically and experimentally in the works of Cliffe et 
a l 111; 13 ■ The transient features of the circular Couette flow were investigated in 
laboratory and numerical experiments (using the commercial program 'Nekton' 
(Fluent GmbH) by Neitzel et al The results of the numerical experiments 
were close to the laboratory experiments. In case of short cylinders with an as- 
pect ratio T ~ 1 (relation of annulus span to the annulus gap width, defined in 
section EJ) usually only one pair of Taylor vortices arises. In general, this case is 
considered to be not really interesting due to the poorness of the patterns to be 
expected. But in the experiments of Aitta et al [7] a series of interesting flow 
patterns for short cylinders was found. A non-equilibrium tricritical points occurs 
when a forward bifurcation becomes a backward bifurcation. In the Taylor - Cou- 
ette problem the tricritical point can be observed for IV = 1.255 in the classical 
case of a pure rotation of the inner cylinder without axial flow as demonstrated 
in [7j and later in The system of two vortices is symmetrical if the rotation 
velocity does not exceed a critical velocity v\. If the velocity of the flow is larger, 
v > vi, then one vortex grows at the expense of the other one. In the case of 
the forward bifurcation the symmetry is broken continuously and in the case of 
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the backward bifurcation - abruptly. In the phase transition language this is a 
bifurcation of the second order or first order correspondingly [9 i . The continuous 
transition was observed for r < IV and the abrupt one for IV < T < Tc, where 
r c = 1.292. 

An exhaustive study of the influence of the finite length effects is done for the 
classical arrangement of two concentric cylinders with a very small aspect ratio 
T G [0.5, 1.6] and a relatively large annulus with a radius ratio (relation between 
the radii of inner and outer cylinders, see section |2J) of i] = 0.667 for moderate 
Reynolds numbers Re G [100, 1500] in the work of Furukava, Watanabe, Toya & 
Nakamura ^0] in the case when the inner cylinder is rotating. The numerical 
investigations were verified by experiments. The three main flow patterns were 
found which correspond to a normal two-cell mode, an anomalous one-cell mode 
and a twin - cell mode. A steady case as well as an unsteady mode of the fully 
developed flow different from wavy Taylor - Couette flow were obtained. 

There exists a series of experimental and theoretical works devoted to the 
linear stability of the Taylor-Couette problem with implosed axial effects. The 
first such analysis was done for the case of axisymmetric disturbances in a narrow 
gap [TT], ^21- In the work of Marques and Lopez the linear stability of the 
flow in the annulus between two infinitely long cylinders, driven by a constant 
rotation and harmonic oscillation in the axial direction of the inner cylinder was 
analyzed using Floquet theory. The axial effects in the Taylor-Couette problem 
were studied in the work [T3| on two different flows. In the first one the axial effect 
is introduced by an inertial axial sliding mechanism between the cylinders and 
in the other one via an imposed axial pressure gradient. In all cases the annular 
gap is much larger as in our case and the radius ratio is 77 = 0.5, . . . , 0.8. The 
stabilizing effect of a periodical axial flow was studied in the works ^H] , ^B] and 
[T7| . In the work of Lueptow [TH] the stability problem of a Taylor-Couette flow 
with axial and radial flow was studied for a relatively large radius ratio rj = 0.83 
and pressure-drive axial flow in the whole annulus. In all these studies it was 
found that the axial flow stabilized the Taylor-Couette flow and enables a lot of 
super critical states with a rich variety of patterns. 

Nevertheless the theory of the Taylor - Couette flow is by far from being 
complete, so for instance the stability problem for the system of Taylor vortices 
and a classical Couette flow by influencing control of the cylinder caps. It is also 
a little bit astonishing that the exact analytical solution for the Couette flow 
profile with boundary conditions on cylinder caps was obtained only recently in 
the work of Wendl [HI]. His results expose a strong influence of the caps on the 
flow profile. The difference between the two profiles with and without finite- 
length effects is approximately a logarithmic function of the aspect ratio F for a 
wide range of radius ratios 77. The solution is given in form of a slowly convergent 
series which contains trigonometric and Bessel functions. To give a qualitative 
analysis of instabilities one takes an exact solution as a base flow to linearize the 
Navier - Stokes equations. This type of solutions is not convenient to pursue the 
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topic further towards analytical studies of the instabilities. 

In the present work we use the theory of small perturbations to study the 
instability of the laminar Couette flow in a short cylinder. The geometry of an 
annular gap and the boundary conditions are given in the section El Our idea was 
to get an analytical formula which is simpler than that in ^Hj but which gives us 
an approximation of the Couette flow with finite-length effects which is as good 
as that found in [TJ] . In the section |3] we describe this solution and discuss the 
region of applicability of the done approximation. 

The solution obtained for a short cylinder is well suitable as a base flow for 
a comprehensive study of instability effects. This study is represented in the 
section The onset of the instability is predicted using the method of local 
analysis, in which the equations governing the flow perturbation are linearized. 
For the very reason that the appearance of an instability depends not only on 
the base flow but also on the amplitude and the type of the flow perturbation 
two different types of perturbations are studied in this section. We consider the 
sensitivity of the Couette flow under small perturbations corresponding to an 
axial flow which is the most important for industrial applications. We take the 
case of a very narrow annulus rj ~ 1 which is typical for bearing or sealing systems 
and in contrary to general expectation we get a rich family of flow states. We 
proved that the axial flow stabilizes the Taylor - Couette flow also in case of very 
short cylinders and narrow annular gap with axial flow directed as a thin layer 
along the surface of the inner rotating cylinder. 

The elaborate series of the numerical experiment was done with help of the 
FLUENT 5.5/6.0. The detailed description of the grids used and others numerical 
facilities is given in the section El 

The change from one laminar flow to another or a transition from a laminar 
flow to turbulence is a function of the base flow considered as well as of the type of 
small perturbations. We found a rich family of typical patterns on large intervals 
of Reynolds numbers. The results of the numerical experiments are discussed in 
the sectional 

2 The main problem 

The study of the Couette - Taylor problem for the parameters defined below was 
stimulated by the industrial project [22] • In many technical systems the problem 
of transmitting energy and momentum is solved using shafts and corresponding 
bushes with or without radial sealing. In a typical situation one is concerned 
with a small Taylor- Couette annulus, narrow slits in between cylinders and with 
axial and tangential oscillations or with a weak axial flow along the surface of the 
inner cylinder. In general, these systems are very complicated and it is necessary 
to study more simple problems in order to isolate a particular mechanism and 
to develop and test hypotheses governing their behavior. We separate here only 
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one of the most typical cases because of its technological importance and, from a 
fundamental point of view, for its rich dynamics and try to describe all possible 
flow states in this case. 

We consider the motion of the fluid in an annulus between two cylinders as 
represented in fig. ^ The annulus dimensions are the gap width d, the annulus 
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Figure 1: The geometry of the annulus. The radii of the inner and outer cylinders 
are equal to ri n ,r out correspondingly, the annulus span is denoted by I, the annular 
gap width is d = r ou t — ri n . The corresponding geometrical parameters of the slits are 
denoted by d and I. Other notations are explained in the text. 

span /, the radius of the inner cylinder r in , the of an outer cylinder r out . The 
corresponding dimensionless geometric parameters are the aspect ratio r = k and 
the radius ratio rj = We use these notations during the whole paper. For 
many sealing or bearing systems T ~ 1 and rj ~ 1. We fixe in all our calculations 
r = 1.02, 77 = 0.95356 . 

We assume that the inner cylinder to turn at the constant rate uj. The outer 
cylinder is at rest, the axial flow has a constant velocity wq and we suppose ri n u ^> 
wq. The equations of fluid motion and boundary conditions will be rendered 
dimensionless using d, as a unit for length. The radius of the inner cylinder will 
be denoted by R = ^ in dimensionless units which will be used in all sections 
of this paper. 

The entrance way and the outgoing way for the axial flow are very thin cylin- 
drical slits with the gap width d and the annulus span of slits denoted by /. The 
slits are long and narrow in comparison to the main annulus i.e. I ~ I, d ^> d. 
The ends of the slits narrow conical to the lower order of scale that means the 
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most narrow part of the slit has the width d, which is very small, d 3> d ^> d. 
A more detailed picture of the slits is given in fig. [7| The studied flow states are 
very sensitive to the kind of perturbations. So we try to construct a situation 
which is most adapted to technical facilities and imitate the sealing arrangements 
with conical tapering slits. 

It is characteristically for quite all technical applications that the sealing 
medium is oil with non vanishing viscosity. We will consider in our work an 
incompressible viscous fluid in the cylindrical annulus. 



3 An approximation of the exact solution in 
case of the rotating inner cylinder 

The classical case of the Couette flow of an incompressible viscous fluid in an 
annulus between two cylinders was investigated analytically in ^H]- in compari- 
son to the famous work of Taylor the author studied cylinders of a finite length 
and got a much more complicated formula which describes the influence of the 
cylinder caps on the flow. As follows from experiments of Kageyama ((201) if the 
Reynolds number grows the influence of the Ekmans layer drops along the cylin- 
der caps. As a result the profile of the solutions differs more and more from the 
pure Couette solution because of the arising of secondary flows in the annulus. 

We will consider the flow in the annular gap with the width d, the span I, and 
radius of inner cylinder r in as on the fig. [T] without slits. We suggest now that the 
caps of the cylinders contact the inner cylinder and that there is no axial flow in 
the annulus. The equations of the fluid motion and the boundary conditions will 
be rendered dimensionless using v = r in uj and P = (v/d) 2 as units for velocity 
and the pressure. The scaling of the velocities and of other quantities will change 
from section to section. 

Due to the axial symmetry two of the components of the velocity vector 
u(r, 9, z) vanish in this case and the velocity vector reads 

u(r,z) = (0,H(r,z),0), (1) 

where Vo(r, z) is the azimuthal component. The cylindrical Navier-Stokes equa- 
tions can be reduced in this model to one equation, 

d 2 V 1 dV ^ Vp 1 d 2 V = 

dr 2 R + r dr (R + r) 2 T 2 dz 2 U 

The boundary conditions can be rendered dimensionless as 

V (r,0) = 0, V(r,l) = 0, V Q (0,z) = l, V Q (l,z) = 0, z e [0, 1], r e [0, 1]. (3) 

We make an additional scaling in z— direction by T to reach the more convenient 
situation that the both variables r and z lie in the interval [0, 1]. The considered 
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case was solved by Wendl jTHj using the Fourier transformation. The solution has 
the following form 

V (r,z)= (4) 

1 (J 1 (/3 m (/? + l))ir 1 (/3 m (i? + r))-ir 1 (/? m (i?+l))/ 1 (/3 m (i? + r)))sin(/j m z) 
4vr ^ (2m - 1) (h(J3 m (R + l^K^R) - K^^R + l))h{(3 m R)) 

with (3 m = ( 2m ~ 1 ) 7r - Here Ji and K\ are modified Bessel functions. This solution 
differs from Taylor-Couette's simple linear solution especially strong in cases of 
short cylinders but remains an exact solution for long cylinders. The series in 
(J3J) converges very slowly and it is difficult to use it for the intended analytical 
studies. This was the reason that we looked for simpler analytical representations 
of this solution or for some good approximated solution. 

We study the case of short cylinders, i.e. T ~ 1, with a narrow gap, r\ ~ 1, 
and R 1. Under these assumptions and taking into account r 6 [0, 1] we can 
simplify the initial equation (J2J) by replacing r^^p, by the constants j^, 
For the simplified equation which has constant coefficients we obtain the solution 
in closed analytical form, 

Vo(r, z) = — arctan ^sinh -1 (p r j — sinh -1 ^— (2 — r )jj sin(7rz)j (5) 

for r G [0, 1] and z G [0, 1]. This is a convenient form of the approximated solution 
for the fundamental mode of the Couette flow profile which we will use in the 
following. 

The agreement between the approximated and the exact solutions is very 
good as can be seen from fig. |2l The approximated solution (jSJ) is stable for small 
Reynolds numbers Re and can be used as a base state in a stability analysis. 

We note an earlier attempt by Vladimirov 21 j to find an approximative solu- 
tion which takes the influence of the caps into account , which, however, differs 
much more from the exact solution than ours. 



4 Stability problem for the Taylor - Couette 
flow under small perturbations 

4.1 Local stability study for a short cylinder 

We consider the axisymmetrical stability problem of the Taylor-Couette flow. 
We assume that the velocity vector u(r, z) which describes the fluid motion in 
the annular gap like on the fig. ^ can be represented as a sum of two vectors 
u(r, z) = (u, v, w) and v(r, z) = (0, V , 0) 

u(r, z ) = v(r, z) + u(r, z) = (0, V , 0) + (u, v, w), (6) 
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a). R = 20, b). R = 10, 

Figure 2: The profiles Vo(r,z) of the exact (dashed line), the approximated 
solutions (continous line) and numerical solution (dash-dotted line) on z— levels 
equal to z = 0.05, 0.25, 0.5 with T = 1 as functions of r. The numerical solution 
was found with program FLUENT 5.5/6.0 in axisymmetric case for Re = 100. 

where Vo describes the fundamental mode of the Couette flow and u(r, z) a pertur- 
bation of this mode. We assume also that the pressure P can be also represented 
as sum P = P + p. 

Under the following assumptions 

I / M I / M |VU| IVvl 

max |v(r, z)\ max |u(r, z)\, max— — — 3>max— — — , (7) 

r,z r,z r,z u r,Z V 

where V denotes the spatial gradient of the corresponding velocity, the system 
of the Navier-Stokes equations can be linearized. The linearized Navier-Stokes 
equations (LNSE) in the axisymmetrical case take the form 

d t u — 2 Re - = — d r p + Am — ^ , (8) 

R + r rF {R + rf V ; 

d t v + Re (d r V u + -^—u + d,V w) = Av - V 

d t w = ~r^d z p + uAw, 

d r u + — h d z w = 0, 

R + r 



where A is the Laplace operator, 
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The Reynolds number is Re = ujri j id and v denotes the viscosity of the fluid. The 
dimensionless variables r, z describing the flow in the cylindrical annulus belong 
to the intervals 

rG [0,1], zE [0,1]. (10) 

Here we use the dimensionless variables as introduced in the section 121 and 
additionally we use the value d 2 / v and Ruo to render time and velocities dimen- 
sionless too. We leave the same notations for the velocities in the dimensionless 
case as above. It means that now |Vo| ~ 1 and from the conditions (|7J) \v\ -C 1 
follows. 

Let's introduce a stream function \1/ as usual in the axisymmetrical case |2II| 

u — t; > w = • (11) 

R + r R+r v ' 

In terms of the stream function \l/ and the azimuthal velocity v the LNSE (|8j) 
turns into a system of two coupled equations 

~ Re(2d z V ^ 2V \ ~ 

dt ^-^\R-rr v+ R-r^ v ) = AAlI/ -(^f' 

dtV + Re((d r V + ^-)4^-d z V ^-) = A-u — 2 • (12) 

VV R + r) R + r z0 R + rJ {R + rf 1 ' 

We remark that the value A\l/ is an azimuthal vorticity Tg defined by 

T s = A\l> = ^d z u — d r w, where A = — — — . (13) 
T 2 R + r 

The basic Taylor- Couette flow is an axial inhomogeneous flow in the case 
of the short cylinder because of the cylinder caps influence. The inhomogeneity 
leads to additional terms which make a difference between the system (JE1) and the 
corresponding linearized Navier-Stokes equations for an infinitely long cylinder. 
Now we will prove that the additional terms causes the propagation of dispersive 
waves in the Taylor - Couette flow. 

We remember that in the case of a short cylinder with very narrow gap we 
have 

T~l, ?7~1, #>1. (14) 

Under these assumptions we can simplify the system of equations (|T2*|) by replac- 
ing the coefficients ^ R ^ n (note r e [0, 1]) by the constants n = 1,2. The oper- 
ators d z and d r can be replaced by the multiplicative operators d r ~ 1/d, d z ~ l/l 
or, in the dimensionless case, by 

dr ' dz T ^ 
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because of the relations (fTUj) and (JHJ). Furthermore we assume 



^ » ^ (16) 
dr R + r K ' 

and neglected -^h in the system We collect now the leading terms of the 

system (I12J1 only and obtain a simplified system of equations 

d t AV -2Red z V v -2ReV d z v = A 2 \I>, (17) 
d t v + 5Red r V d z ^ -5Red z V d r y = Av, 

where no-slip boundary conditions u\s = v\g = w\s = on the inner surface S 
of the annulus are adopted. In the system (JT7J) we denoted A = <9 2 r + fsd% z , and 
the constant 8 = 4 is small, 5 << 1. 

The standard, straightforward way to investigate the stability problem would 
consist in finding the exact solution of the eigenvalue problem corresponding to 
the system of equations (|17)1. Then the perturbations u(r, z) can be represented 
as a series over the eigenfunctions. Alternatively, the system could be solved nu- 
merically whereas a spatial discretization of the problem may be accomplished by 
a solenoidal Galerkin scheme. In our case these procedures may be accomplished 
numerically only and this is outside the scope of the present work. 

Nevertheless we can give good qualitative characteristics of the flow in the 
annulus if we use a local stability analysis. This method was first developed in 
the analysis of an inhomogeneous complex astrophysical flow [211,125] • It give 
us the possibility to establish the main physical processes of the considered fluid 
flow system [2H], [23- 



The main idea of the local stability analysis method also named a method 
of infinitesimal perturbations can be described as follows. We suppose that the 
amplitude and the gradients of the main state (the fundamental mode of the 
Couette flow in our case) is preserved under the perturbations. This means we 
keep the values Vo, d z Vo and d r Vo unaltered or assume that they are constant. For 
all components of the perturbation vector u(r, z) = (u, v, w) and their derivatives 
Taylor series are used which are truncated after the first two terms. From the 
physical point of view we assume that the perturbation scale is much smaller 
than the typical scale of the axial inhomogeneity in the Taylor- Couette system. 
This is the locality condition. As usual we represent the perturbation as a plain 
wave, 

~ \& exp (yt + % ik r r + k z z)) , v ~ v exp (ft + i (k r r + k z z)) , (18) 

where ^ and v are the amplitudes of the waves. The vector of wave numbers 
k = (k r , k z ) defines the spatial direction of the plane wave propagation. The 
locality condition requires that k r ,k z ^> 1. In ()18j) . 7 has the dimension of a 
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frequency and in the stability analysis it is called the time increment (or simply 
increment). In the case of an intrinsic instability of the investigated system 7 
will be positive. This means that the amplitude of the plane wave will infinitely 
increase with time. This is the unstable case. 

The substitution of the plane wave (fTHj) into the linearized system (fTTj) results 
in a system of algebraic equations for the increment 7, 

-v Re 

- A; 2 (7 + - 2— (d z V + iV k 2 )v = 0, 

(j + k 2 )v + i6Re(d r Vok z -d z V kr)$ = 0, 

where k 2 = k 2 + T~ 2 k 2 z . This system can be solved explicitly and we obtain two 
solutions for the increment 7, 

2 Re(25) 1/2 \k z d r V - k r d z V \ 1/2 , 2 2 T/ x2\i/4 
7± = -k ± ^ [k;V + (d z V ) ) x 

(S)—e— (©))■ 

As a consequence, in the system two different regimes of motion can exist. 

The motion in the annulus is strong inhomogeneous in z— direction because 
of the influence of the cylinder caps. Let us now look for local increments on the 
middle surface r = 0.5 on three different lines z = 0.25, 0.5, 0.75 . For z = 0.5 the 
increment does not have any imaginary part and on this line the system does not 
show any oscillations. In this case the real part of the increment is represented 
in fig. El On the lines z = 0.25, 0.75 the increment has an imaginary part and 
the system has oscillations with frequencies which are directly proportional to 
the gradient of the undisturbed state. 

For d z Vo = we have two real solutions for the increment and, consequently, 
two branches of non-oscillating solutions. 

The two branches of the increment given by Eq. (fTU|) may have also an oscillating 
character. For d z V$ 7^ we deal with two oscillating solutions corresponding to 
waves of the inertial type. If additionally k z d r Vo — k r d z Vo = holds than these 
two oscillating branches have an intersection point. In the intersection point 
7± = —A; 2 we have to do with a time decreasing solution. On the other hand the 
perturbations with k z d r V — k r d z V ~ have the possibility to exchange energy. 



4.2 Local stability of the Couette flow under a small axial 
perturbation flow 

The study of the stability problem of the Couette flow was done up to now for 
arbitrary small perturbations. Let's consider the effects of an axial flow near the 
surface of the inner cylinder. This flow can be generated by the pair of Taylor 
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Figure 4: The imaginary part Imj(k r ) 
of the increment 7 as a function of k r 
for z = 0.25, r = 0.50, k z = -2.00, 
Re = 110. 
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Figure 5: The real part Rej(k r ) of the 
increment 7 as a function of k r for z = 
0.50, r = 0.5, fc, = -2.0, Re = 110. 



Figure 3: The real part Re~j(k r ) of the 
increment 7 as a function of k r for z = 
0.25, r = 0.50, k z = -2.00, Re = 110. 



vortices or an axial flow through the narrow slits in the caps of the cylinders near 
the surface of the inner cylinder as shown in figure ^ 

In order to simplify the study we suggest that the axial velocity w(r,z) 
has a nonzero value w(R, z) on the surface of the inner cylinder. Addition- 
ally we suggest that the boundary condition is homogeneous in z— direction, 
i.e. w(R, z) = ( = const holds. Here we can neglect the dependence on the 
z— coordinate. The constant constituent of the axial velocity distribution can be 
removed by an Galilean transformation z' — z — £ (1 + R) t; z' — > z. The bound- 
ary condition w(R,z) = ( ^ provides in the narrow gap approximation the 
linear velocity profile in the r— direction 

wK£(l-(r-R))=t(l + R)- (r. (20) 
Now we can rewrite the LNSE (JHJ) in the form 

a Air, a a,t, ■ CRec^E 2Red z V 2ReV ~ A* 

d t A^> - (Rerd z A^ + - 9 rV - — rd z v = AAW 



(i? + r ) 2 r 2 (i? + r) r 2 (i? + r) (i? + r) s 

3 t v - C Re rd z v + Re (d r V& + t^-^ - Re d z V -^- = Aw 
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We will simplify this system of equations. We remind that we investigate the influence 
of small perturbations and from the conditions QC^CVoor^^Cl follows immediately. 
We expand the coefficients in the system into a series with respect to the parameter 
5 = 1/R and keep the main terms only. The approximated system for the stream 
function and for the azimuthal velocity v has the structure 

2Re 2Re 

d t Ay-(Rerd z Ay + (6Red z y--^-d z V v--^rV d z v = A 2 *, 

d t v-CRerd z v + SRed r V d z -$ -5Red z V d r y = Av. (21) 

As in the general case in the previous point H~T1 we will study the stability problem 
under small local perturbations. So the perturbation scale is much smaller than a 
typical scale of the inhomogeneity in the system. This locality condition implies to the 
inequalities 

k r > l,k z > 1. 

In the same way we assume that the values Vo,d z Vo and d r Vo change very slowly 
in comparison to the perturbations so that they can be represented by Taylor series 
truncated after the first two terms. As usual we assume that the perturbation u(r, z) = 
(u, v, w) has the form of a wave with time dependent amplitudes and with a time 
dependent wave number in r— direction, 

= * (t) exp (i {k r (t) r + k z z)) , 
v = v(t) exp (i (k T (t) r + k z z)) , (22) 

where \I/ is the stream function IjllJI and v is the azimuthal velocity. The substitution 
of the formulae (|22|) into the system of equations ()21j) leads to the simplified system of 
equations 

-8 t ((A; 2 (t) + k 2 jT 2 ) $ (i)) - i (d t k r (t) -(Rek z ) r (k 2 r (t) + k 2 z /T 2 ) $ (t) + 

iC8Rek z * (t) - ^ (d z V + tV k g )v(t) = (k 2 (t) + k 2 /T 2 ) 2 * (t) , 

d t v (t) + i (d t k r (t) - C Re k z ) rv (t) + i6 Re (d r V k z - d z V Q k r ) $ (t) = 

-(k 2 (t) + k 2 /F 2 )v(t). (23) 

The solution of this system can be found for all r if the condition 

d t k r (t)-(Rek z = Q, (24) 

is fulfilled because only two terms in the system (|23[) have coefficients which depend 
on the variable r. The equation (|24jl can be easily solved and we obtain 

kr(t) = k r (0) + (Retk z . (25) 

The radial component k r of a wave number vector k increases or decreases linearly 
in dependence of the sign of the axial component k z of the wave number vector. In 
another words we can explain the action of the axial shear as a change of the radial 
scale. 
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The time dependence of the wave number k r {t) leads to an algebraic non-exponential 
evolution which can be described by the system of equations 

-8 t ((fc 2 (t) + k 2 z /T 2 ) $ (i)) + i(5Rek z $ (t) - ^ (d z V + iVb**)t> (*) = 

(^(tj + ^/r 2 ) 2 ^), 

d t v(t) + i5Re (d r V k z - d z V k r ) § (i) = - (A; 2 (t) + fc 2 /r 2 ) u (t) . (26) 

Perturbations as studied here were considered in fluid dynamics since Lord Kelvin, but 
it has been realized only recently that these perturbations can play a trigger role in 
the subcritical turbulence transition. These untypically properties of the flow resulting 
from the fact that the corresponding LNSE operator is non-self-adjoint for shear flows 
such as a plane Couette flow [2S] , j^H] , a Poiseuille flow jSOj , a counter-rotating Taylor- 
Couette flow [HJ, jS2], a boundary layer [HH], a Kepler rotation [2B]. Apart from 
the modes of the discrete spectrum there exist plenty of modes of the discontinuous 
spectrum. The algebraic time evolution of these modes can achieve large amplitudes 
which can reach hundredfold of the initial values of amplitudes. Another remarkable 
property of these modes is their trait to redistribute energy and momentum between 
different types of perturbations [21], 

To give the qualitative description of an impact from an axial flow on the stability 
of the Couette flow we will find the increment for a fixed instant time. Analogously to 
((TH|) we have the representation 

~ exp {"ft) , v {t) <~ exp {"ft) (27) 

for the time dependent amplitudes of the wave (|22|). If we virtually hold unaltered 
the time-dependence of k r {t) we can consider the stability analysis local in time. The 
substitution of the formula (f2*7|) into the system of equations (|2f)|) leads to an algebraic 
equation for the increment. We get now two solutions of this system for the increment 

7 

_ 2 i(5Rek z 
7± " k + 2k? ± 



Re^ 1 / 2 (64/fc 4 A 2 {d z V ) 2 + {( 2 5k 2 z + 8k 2 k z V A) 2 ^ 
A = k z d T Vo — k r d z Vo, <3? = - arctan ' " 



exp H$) , (28) 



2 \( 2 5k 2 z +8k 2 k z V A 

For small (Re the dispersion curves have the same form as considered in the previous 
section (without any axial flow). In the intersection area k z d r Vo — k r d z Vo ~ we have 
now 

2 i(5Re k z 2 / on x 

7+ = -k H ^2 — , 7- = ~k ■ (29) 

This dispersion equation is similar to the equation for the Rossby waves for inhomoge- 
neous rotation of the Earth or othetrs systems. The wave number k r {t) is no longer fixed 
as in the previous case Q18|) but it has a time dependent value. As a consequence, the 
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point Irwy (k r (t)) or Re^f (k r (t)) will move along the branches (similar to the branches 
on the fig. El - fig- EJ) with time. Chagelishvili and Chkhetiani |34j have proved that 
under the influence of the shear the flow frequencies of the Rossby waves and inertial 
waves will be modified so that they take the same values. This circumstance leads to 
the energy exchange between different modes in the flow. The same mechanism of the 
energy exchange is characteristic for the system (j26|) . 

Suppose that we start the motion near the intersection point. One mode is unstable 
and starts to grow while the second mode will dissipate. In the intersection point their 
frequencies are very close to each other and a resonance interaction may take place. In 
this case the dissipative mode loses the main part of its energy. Thus, we have found 
a mechanism for the instability which may be responsible for transport and exchange 
of energy and angular momentum. 

If we analyze the solutions of the system (|26j) we notice that axial flow leads to a 
strong differentiation in the energy growth of the perturbations in dependence on the 
axial level. This dependence is represented for two Reynolds numbers in fig. El This 
asymmetry can be the seed for the symmetry breaking in the Taylor-Couette system 
in short cylinders. 



Figure 6: The time dependent energy growth for T = 1, ( = 0.1 , r = 0.5 , k r (0) = 1.0, 
k z = -1, $(0) = 1.0, u(0) = 1 for different levels of *, z = 0.25 , 0.50 , 0.75. 

We can assume that the same mechanism work also without any axial flow in the 
system. Consider two symmetrical counter rotating Taylor vortices in the annulus. 
Then we have in the annular gap a radial distribution of the axial velocity component 
which is oppositely directed in different parts of the annular gap. If now a small 
perturbation amplified the circulation in one of the vortices for a small interval of time 
then radial gradients of the axial velocity will show up in different parts of the gap. 
Therefore the degree of growth of perturbations in the continuous spectrum will be 
different in both parts of the gap and the faster growing perturbation can trigger an 
irreversible process of symmetry breaking. 





b) Re = 33 
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5 A symmetry breaking and helicity conserva- 
tion law 



The classical Taylor - Couette problem with pure rotation of one or both cylinders 
posses Z 2 and SO(3) symmetries. The implementation of an axial flow immedi- 
ately destroyes the Z 2 mirror symmetry. The axial perturbation flow amplifies in- 
homogeneously different intrinsic modes in the studied states as we proved in the 
previous section. It results in a reach variety of super critical stationary and quasi 
periodical states. This symmetry breaking for short cylinders can be considered as a 
generation of certain helicity states. 

The helicity is a topological characteristic of the flow and it is the second inviscid 
integral invariant for the Euler equations. The helicity is defined by 

H = J u • [V x u] dV, (30) 
v 

where the vector T = V x u defines the vorticity of the flow. The helicity is a pseudo 
scalar because of 

W(-r) = -W(r). 

It is well known that the helicity is an important quantity in the study of dynamics 
and stability of complicated vortex structures and flows. The states with maximum of 
helicity correspond to the states with minimum of energy |36| . |37j . |38j . 

In the same way as we get the evolution equation for the energy of the system 



d t £ = jp- u n dS -v J [V x u] 2 dV 



S V 
we obtain the following evolution equation for the helicity, 



d t H = j> [V x u] n f y - dS - v J [V x u] • [V x [V x u]] dV, 

S ' V 



(31) 



where n is the unit normal vector orthogonal to the boundary surface S. 

For the inviscid fluids the right hand side of equation (|31j) reduces to the first term. 
The helicity conserves if the normal component of the vorticity on the boundary surface 
disappears, [V x u] n = 0, or if we have a potential flow. Let S include a bounded 
volume of an ideal fluid with T / and TL ^ then the helicity will be conserved 
by the motion. For very small viscosity v ^ the helicity is no longer conserved but 
it changes very slowly because of the negligibleness of the second term in H31|) . As a 
consequence, for large Reynolds numbers we can treat the helicity as an approximately 
conserved quantity. 

In a classical Taylor-Couette system with mirror symmetric boundary conditions we 
expect an even number of counter-rotating vortices. The growth of anomalous modes 
and the phenomenon of appearing of time-dependent states can be considered as a 
mirror symmetry breaking process and a generation of states with a definite sign of the 
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helicity. Let us represent the helicity in the cylindrical coordinate system and reduce 
the formula to the case of short cylinders. The components of the vorticity T are 

u 

T r = —d z v, Tg = d z u — d r w, T z = d r v H — . (32) 

r 



Using the expressions (|32j) we obtain for the helicity 7i in the axisymmetrical case the 
convenient representation 

TL = 47r / / (vd z u + wd r v H ) rdrdz. 

r J 

For a narrow annular gap with rj ~ 1 and ]? > 1 we can approximate the formula for 
the helicity by 

U 2* 4ttR J j (yT e + drdz ra 4ttR J j vT drdz. 

s s 

The axial flow generates a constant stream of helicity and angular momentum through 
the annulus. It is an important source of the symmetry breaking in the system. 



6 Numerical facilities 

Major advances in computer technologies and numerical techniques have made possible 
to propose an alternative or at least a complementary approach to the classical and 
analytical techniques used in laboratory experiments. Computational Fluid Dynamics 
becomes part of the design process. But discontinuities or deep gradients lead to 
computational difficulties in the classical finite difference methods, although for the 
finite element methods often a suitable variational principle is not given. Therefore, 
we consider the finite volume methods, which are based on the integral form over a cell 
instead of the differential equation. Rather than a pointwise approximation at the grid 
points, we break the domain into grid cells and approximate the total integral over 
each grid cell. It is actually the cell average, which equals this integral divided by the 
volume of the cell. 

We consider the axial symmetric case as well as the full 3D flow in the annulus. 
The numerical method allows to choose a so-called segregated or a coupled solver. The 
simulation process consists of (see e.g. |3*H] ) 

• division of the domain into discrete control volumes (construction of the compu- 
tational grid), 

• integration of the governing equations (in axial-symmetric or in full 3D form) on 
the individual control volumes to construct algebraic equations for the discrete 
variables (velocity, pressure, temperature and other parameters), 

• linearization and solution of the algebraic system to yield updated values of the 
variables. 
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Both the numerical methods of the segregated and the coupled solvers employ a similar 
finite volume process, but the approach used to linearize and solve the equations is 
different. 

The FLUENT 5.5/6.0 [101 simulation and post-processing program was applied for 
the numerical calculation of the flow. The unsteady 2D axisymmetric or 3D version 
of the Navier-Stokes equations is solved as described using finite volume discretization 
on a structured quadrilateral grid (fig. EJ). The finite volume method has a first order 




Figure 7: 2D structural grid for the annulus and the slits with 7690 quadrilateral 
cells. 

upwind scheme and constant relaxation factors, while the SIMPLE procedure for the 
pressure-velocity coupling is applied. 

The grid is generated with the help of GAMBIT program and characterized by a 
large difference (100 times) between boundary edges. But it was possible to save the 
grid structural and quadrilateral. The smallest inflow and outflow faces of the grid 
are divided uniformly by 20 nodes, whereas the nodes in the bigger geometry of the 
axial flow are a more concentrated near the boundary because some vortices can arise 
there if the Reynolds number is large |3J|. The major volume of the annular gap is 
divided 50x70 nodes, uniformly in the axial direction and linear condensing in the radial 
direction. 

As next, we set up a 2D axisymmetric model, using a rotating reference frame. The 
time step is selected to be ^ of a single rotation. It helps to determine the flow motion 
state in the quarters of the rotation circle {0 = 0,7r/2, tt, 3tt/2). Each time step includes 
30 iterations or fewer if the solution converges with the residual of order O(10 -5 ). 

The plane of the 3D grid that corresponds to 2D axisymmetric model have to be 
coarser, because of the large diameter of the inner cylinder and therefore of the large 
number of cells in the projection on the (r, #)-plane. These are only 10 nodes on the 
inflow and outflow faces in the conical tapering of the slits and 25x30 in the annulus of 
the slits (fig. El). It results in 1087040 hexahedral cells, that is close to the limit for a 
computation and post-processing on a single processor in an acceptable time. Another, 
more complicated 3D simulation would be possible using parallel computing |32]. Such 
a program is used to being written with the Help of Message Passing Interface (MPI) 
|43j or some other specific approaches [H]. But this is the aim of further investigations 
and developments. 
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Figure 8: 3D grid with 1087040 hexahedral cells. 



Some test simulations showed, that the 3D case is not sensitive enough to determine 
the very small scale instabilities correctly. Although, the results of the 3D and 2D 
axisymmetric simulations on the same grid are equivalent. It was verified using the 
same coarser grid for the 2D case or by solution of the stationary problem without 
instabilities. A finer grid for the 2D axisymmetric case, by contrast, does not give more 
substantial changes to the behavior of the unstable flow. So, the 2D axisymmetric case 
(fig. EJ) with 7690 quadrilateral grid cells is chosen and can be considered without loss 
of accuracy for all computations in that geometry. 

7 Pattern structure in the annulus under influ- 
ence of the axial flow 

In numerical experiments which simulated the flow in the annular gap as in fig. ^ 
and more detailed in fig. both axial symmetrical and 3D cases were studied. 3D 
calculations had been done prepared for a large number of examples. But it became 
apparent that no torsion or other remarkable 3D effects arose in the system for Re < 
6000. Therefor we used for further numerical investigations axisymmetrical case which 
saved a lot of time and computer capacity. As a graphical representation of the results 
in 3-D cases on a paper sheet does not show the interesting region in sufficient detail so 
that we represent only one fig. El which gives an impression, that the processes under 
investigation are indeed mostly axisymmetrical. 
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Figure 9: 3D simulation of the two symmetrical stationary vortices. Only one 
of the vortices is represented to avoid an overload of the left figure. The counter 
lines of the velocity magnitude are given in 3D case on the left and in cross 
sections on the right right by 6 = 0,7r/2,7r, where v = 0.1(10~ 3 — ), w = 0.016 , 
Umax = 0.060 , Re = 1256.0 , Re ax = 0.1 . 

In all investigations the inner cylinder rotates with the same azimuthal velocity 
Vq = 2nu;ri n on the surface of the inner cylinder. We used this value to render other 
velocities dimensionless. The weak axial flow is directed on all figures from left to right. 
The entrance and outgoing slits in the ends of conical tapering have the radius ratio 
fj = 0.99992. 

We found that there exist four different flow states with typical patterns in depen- 
dence on the relation between azimuthal and axial Reynolds numbers. The azimuthal 
Reynolds number Re is defined by Re = Vq d/v and the axial Reynolds number by 

R^ax — tVmax d/V. 

We studied a pure rotation and then switched on an axial flow to investigate the 
sensitivity of the flow patterns to small axial perturbations. The pattern structure 
is very sensitive to the appearance of an axial flow. We looked for very weak axial 
flow with Re ax <C Re nevertheless the pattern structure changed immediately under 
the influence of the axial flow. We studied the states distribution in the regions Re £ 
[0,12000] and Re ax £ [0,10]. The results for the most interesting part are that for 
Re G [0, 2500] and Re ax € [0, 1] and they are represented on the schematic diagram 

fig-EB 

Let us consider the pure rotation case without any axial flow. The most probable 
state for a short cylinder with an aspect ratio T ~ 1 can be described as follows. The 
nearly quadratic cross section of the annulus is filled with two counter rotating stable 
Taylor vortices as represented on the fig. 111! We proved that this structure remains 
mirror symmetric for a large region of azimuthal Reynolds numbers, Re G [0,2500]. 
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Figure 10: The pattern distribution in dependence on the azimuthal Re and 
axial Re ax Reynolds numbers. Two Taylors vortices in the annulus are denoted 
by □, one stable vortex by •, natural oscillation states of two vortices by 
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Figure 11: Two symmetrical stationary Taylor vortices which are denoted by □ 



on the figUUl Here is v = 0.8(10" 4 ^), w = 0, u max = 0.07 , Re = 1570, Re ax = 
0. 



The investigated states are denoted in fig llUl bv boxes, □, on the horizontal axes. All 
of them have quite the same structure - two stationary Taylor vortices in the annulus 
up to very high azimuthal Reynolds numbers. 

An enlargement of the azimuthal Reynolds numbers leads to the second state in 
the annulus. For very large azimuthal Reynolds numbers Re > 4000 the intrinsic 
instabilities give rise to perturbed states with a lot of small vortices. This flow pattern 
can be interpreted as a region with a fine-grained structure overlaid with a large scale 
structure - three Taylor- like vortices (see fig. I12[) . Consider now the onset on an axial 



Figure 12: A lot of fine-grained structures of small vortices overlaid by large scale 



vortices. Here v = 0.1(10~ 4 ^), w = , u max = 0.04 , Re = 12560 , Re ax = 0. 



flow with relatively small velocity. The existence of the axial flow can be considered 
as a small perturbation of the base flow with pure azimuthal velocity. In the region of 
the azimuthal Reynolds numbers Re < 1500 the axial flow promoted one of the counter 
rotating Taylor vortices in appropriate direction. As result the annulus will be filled 
also for very small axial amplitudes of the axial flow with one stable vortex (see fig.H3|l. 
We see the abrupt change of the patterns structure from two stable Taylor vortices to 
one by switching on an axial flow. This is the third pattern structure observed in this 
system. It is denoted by • on the figlTUl 

In the region on the Reynolds numbers He az ~ 1500, Re ax ~ 0, 4 there exists a 
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branch point in which three possible flow states meet. For all axial Reynolds numbers 
0.05 < Re ax < 4 and Re < 1500 it is the first branch for which is characterized by one 
stable vortex in the annulus (fig|13|). For 1250 < Re < 2400 and for the small axial 




Figure 13: One large stationary vortex filling out nearly the whole annulus 
(denoted on the figHOby •). Here v = 0, 1(10" 3 ^); w = 0,05 ; u max = 0,04 ; 
Re = 1256 ; Re ax = 0, 32. 



Reynolds numbers < Re ax < 0.64 a second branch with two stable Taylor vortices 
appears (fig llll) . For Re Qv2 > 1500 and 0.3 < Re ax < 0.9 the third branch appears with 
two quasi periodically natural oscillating vortices (fig lTH) . 




a) t\ = 0.402 b) ti = 0.427 (minimum energy) 




c) 13 OA is ^ £ 4 _ o 4gg (maximum energy) 

Figure 14: The natural quasi periodic oscillations of a two vortex state in the 
annulus (denoted on the fig. El by -fc) in four different moments of dimensionless 
time. Here is v = 0.3(10" 4 ^), w = 0.30 , u max = 0.05 , Re = 5108 , Re ax = 6.333. 

This branch representes non stationary states. The energy and momentum of the 
system show also typical quasi periodical oscillations. The time dependence of the 
energy is given on the fig. 1151 The time dependence of the angular momentum is of the 



23 



same type. We rendered the energy dimensionless using its value for small Reynolds 
numbers, Re = 251.2, Re ax = 0.02 . In fig. ED we also see that maximal values for 
the energy are smaller in the perturbation state as in the non perturbed flow. We will 
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Figure 15: The time dependence of the total energy E(t) for the flow by Re = 
5108, Re ax = 6, 333. 

distinguish two different types of energy in the system - the rotation or toroidal energy, 
defined by 

Etor- = 2vr J v 2 rdrdz 
and the energy caused by axial and local radial flow in the system 

Epoi = 27r J (u 2 + w 2 ) rdrdz, 

named also poloidal energy. Correspondingly with Reynolds numbers the toroidal 
energy is three order larger as the last one. But the intrinsic instabilities in the system 
lead to steady energy exchange between different motion types. This exchange we 
demonstrate by a parametric representation in the plane of both energies in fig. Ilfil 

The flow field between differently shaped bodies of revolution were studied exper- 
imentally or theoretically rarely. The experimental work devoted to the case which is 
somewhat analogously to the motion in the conical slits at the caps of cylinders in our 
case was done by Wimmer |45j . For a rotating cylinder in a stationary cone he proved 
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Figure 16: The parametric representation of the energy exchange between E tor 
and E po i in the plane of both energies by Re = 5108 , Re ax = 6, 333. 

that due to the non constant gap different types of motion can coexist in the gap. By 
the growth of the Reynolds numbers the first vortices appear at the location of the 
largest gap size and are deformed, having larger axial extension. The numerical studies 
in our geometry confirm this statement as can be seen on the fig. 1121 



8 Conclusions 

The convenient asymptotic representation for the exact steady solution of the Taylor- 
Couette problem was found in the case of closed short cylinders T ~ 1 (section EJ) . The 
found approximated solution © will be an exact solution for infinitely large viscosity 
or for Re — ► 0. For finite Reynolds numbers a flow in a thin layer along the caps of 
the cylinder can arise. We estimated the radial and axial components of its velocity. 
At first blush in the case of the strong nonlinear motion the estimation u/Vq ~ R^ 1 ^ 2 , 
w/Vq ~ -R -1 / 2 must hold. This follows from the comparison of the nonlinear terms in 
the Navier-Stokes equations. But the numerical computations for the short cylinders 
with R S> 1 and large Reynolds numbers lead to another estimation, u/Vq ~ R~ l , 
w/Vq ~ R , which can be obtained from the Navier-Stokes equations in the case of 
very large viscosity. In fig. El we represent the radial component of the dimensionless 
velocity only. For the other component the pattern distribution is quite the same. 
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Figure 17: The distribution of the radial component of the dimensionless velocity 
max u in the Taylor-Coutte flow for different states in dependence on the Reynolds 
numbers Re and Re ax . Two Taylors vortices in the annulus are denoted by □, 
one stable vortex by •, natural oscillation state of two vortices by 
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From the estimations it follows that the azimuthal component of the velocity has 
the main influence on the stability of the system. Consequently, we can use for our 
qualitative stability analysis the found approximated solution. 

The stability analysis done in section 0] showed that the inhomogeneity of the base 
flow in the radial direction as well as the axial flow cause a perturbation of the Rossby 
waves type. We proved that the Z 2 — symmetry breaking due to the axial flow in the 
Taylor-Couette flow is in the same time the helicity symmetry breaking (section . 
The states with non-zero helicity can accrue in the system (section Q) . The energy 
oscillations in the case of the unsteady pattern with natural oscillations as in fig. ED 
correspond to the helicity oscillations in the system. This can be seen from fig. 1141 

The weak axial flow in an annulus between the rotating inner cylinder and the fixed 
outer cylinder has several important engineering applications in particular if the axial 
flow is directed as a thin layer along the surface of the rotating cylinder. We studied 
this case and proved that the axial flow stabilizes the motion and that in a system with 
very narrow annulus different super critical states are possible (section EJ). We remark 
that experiments in short cylinders with an aspect ratio of T ~ 1 without any axial 
flow showed a poor variety of states 0Z1- But for smaller aspect ratio all the states 
which we described in our case can be observed. This behavior of the Taylor-Couette 
system can be interpreted as follows. The influence of the weak axial flow along the 
surface of the inner cylinder can be considered as a change in the geometry, i.e., as 
a shortening of the axial length of the cylinder. This interpretation is supported by 
simple estimations. 

For technical applications it is important to describe not only the velocity field but 
also the pressure distribution in the device. The pressure differences in the annulus 
coming from the vortex motion can be estimated as R~ l of the pressure difference 
corresponding to the axial flow. 
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